{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "99bbaa7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "from sympy import symbols, cos, sin\n",
    "from sympy.plotting import plot_parametric\n",
    "from math import sin,cos\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73e974e3",
   "metadata": {},
   "source": [
    "# Лабораторная работа 1,Задание 1\n",
    "Нахождение с помощью расширенного фильтра Калмана(EKF) отслеживаие ориентации и положения робота по наблюдениям ( подбробности модели в условиях лабораторной работы lab1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "62668378",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Радиусы колес\n",
    "B = 4\n",
    "# Шаг дискретизации\n",
    "T = 0.003\n",
    "# Угловые скорости колес\n",
    "W_L = 5\n",
    "W_R = 3\n",
    "#Коэфиценты для угловых скоростей\n",
    "a = 0.4\n",
    "c = 0.4\n",
    "# Угловые скорости колес\n",
    "def u_R(k): return a * (k*T);\n",
    "def u_L(k): return c * (k*T);\n",
    "# Линейные скорости колес\n",
    "def S_L(k): return W_L * u_L(k);\n",
    "def S_R(k): return W_R * u_R(k);\n",
    "# Скорость передвижения\n",
    "def s_t(k): return (S_L(k) + S_R(k))/2;\n",
    "# Скорость вращения\n",
    "def s_r(k): return (S_R(k) - S_L(k))/(2*B)\n",
    "# Функции движения\n",
    "def r_t(t): return (-0.5*  W_L *c* t * t+ 0.5 *W_R * a *t * t)/(2 *B)\n",
    "def x_t(t): return (W_L*c + W_R*a)*B*sin(((W_L*c - W_R*a)*t*t)/(4*B))/(W_L*c - W_R*a)\n",
    "def y_t(t): return -(-W_L*c - W_R*a)*cos(((W_L*c - W_R*a)*t*t)/(4*B))*B/(W_L*c - W_R*a)\n",
    "#Среднеквадратические отклонения шумов показаний\n",
    "sigma_r = 0.62\n",
    "sigma_x = 0.88\n",
    "sigma_y = 1.09"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dae3e3a5",
   "metadata": {},
   "source": [
    "Генерация истинной траектории дивжения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "dd291466",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb710dd8250>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhzUlEQVR4nO3deXhUZZ728e8vGyQBAiRhJwk7KAqYEPYAbqDtQrsjIAqIGy5o29OvPdPtdDs902Pb7oKAiAqCoo3iDtKyrwEB2fcAEQj7voTkef9InGEwkZCq1KlK7s915SJ1UpVzeyR3Dk895zzmnENEREJPmNcBRESkdFTgIiIhSgUuIhKiVOAiIiFKBS4iEqIiArmzhIQEl5KSEshdioiEvKVLl+5zziWevz2gBZ6SkkJmZmYgdykiEvLMLKuo7RpCEREJUSpwEZEQpQIXEQlRKnARkRClAhcRCVEqcBGREKUCFxEJUSFR4LM27GXs3K0cPpnrdRQRkaAREgX+7Zo9/OnzNXT8ywx+9/FKVmUf9jqSiIjnLJALOqSlpbnSXom5Kvsw4xdm8cnybE7l5tO2YXX6d0zmhsvrUjky3M9JRUSCh5ktdc6l/Wz7hQrczMYCNwA5zrnWhdtuB54FWgHpzrkStbIvBf6Twydz+ceynby3MIste49TPSaS21Mb0K9DMikJsT59bxGRYORLgWcAx4B3zynwVkA+8Cbwm0AW+E+ccyzYsp/xC7OYtnoPZ/Md3ZolMKBjMle2rEVEeEiMDomIXFBxBX7Bm1k552abWcp529YWflO/BbxYZkbnJgl0bpLAniOnmLR4BxMXb2foe0upF1eZvulJ3JnekFpVK3uWUUSkLJVoDLywwD//6Qz8nO0zucAZuJkNBYYCJCUlpWZlFXlTLb84m5fPt2tzmLAoizkb9xERZvRuXYeBnVNIS67h6S8cEZHSKvUZuK+cc6OAUVAwhFKW+4oID6N36zr0bl2HLXuPMWHRdiZn7uDzlbu4pG417u2Swk1t6ulNTxEpF8rtQHHjxCr82w2XsPCZq/jLry/jbH4+v/1oJZ3+cwZ//Xod2YdOeh1RRMQnAV3QwQsxURHc3SGJvukNWbBlP+/M38abszbz5qzN9Lq0YHilQ6OaGl4RkZBzwQI3s4lADyDBzHYCfwQOAK8CicAXZrbcOderLIP66tw3PXcePMH4hduZtGQ7X63aTcs6VRnYOYU+besTHaXhFREJDSFzIU9ZOJWbx6fLsxk3P4u1u44QFx3JXe0b0r9jMg1rxngdT0QE8GEeuD8FW4H/xDnHkm0HGTd/K9+s3oNzjqta1ebezil0bhKv4RUR8ZRns1BCgZmR3qgm6Y1q8uOhk0xYlMXExTuYvmYPzWtXYVCXRvRpV1+zV0QkqOgMvBincvP4fOUuxs7dyppdR6gZG0X/Dkn075Ssi4NEJKA0hFJKzjkWbjnAW3O3MmPdHiLDwrixTT0Gd23EJfWqeR1PRCoADaGUkpnRqUk8nZrEs3Xfcd6et5XJmTv5eNlOOjeJZ3DXRvRsUYuwMI2Ti0hg6Qy8FA6fyGXiku28M38buw6fonFCLPd1SeHW1AbEROl3ooj4l4ZQykBuXj5frdrNW3O3smLHIeKiI+mbnsTAzsnUjYv2Op6IlBMq8DLknGPZ9oO8NXcrX6/aTZgZ119Wl8FdG9GmYXWv44lIiNMYeBkyM1KTa5KaXJMdB04wbv42Pliyg6krfiS9UU0eyGiscXIR8TudgZeRo6dy+WDJDsbO3cqPh0/RrFYV7s9ozM1t61EpQvPJRaTkNITikdy8fL5YuYuRszazbvdRalWtxH1dGnF3hyTioiO9jiciIUAF7jHnHHM27mPU7C3M3bSPKpUi6JvekEFdG+kNTxH5RSrwILIq+zCjZm/hix92YcBNbepxf0ZjWtXVhUEi8nMq8CC048AJxs7bygdLdnDiTB4ZzRN5MKMxnXQDLRE5hwo8iB06cYbxC7MYNz+LfcdO07p+NYZmNOH61nWICC+3iyaJSAmpwEPAqdw8pnyfzejZW9iy7zhJNWN4sHsTbk2tr5krIhWYCjyE5Oc7pq3Zw4iZm1ix8zC1qlZiSLdG3N0hmSqVNHVfpKJRgYcg5xzzN+/n9e82MX/zfuKiIxnYOYX7OqdQIzbK63giEiAq8BC3fMch3vhuE9PW7CEmKpy+6Unc360xdeJ0b3KR8k4FXk5s2HOUkTM38+mKHwkzuPWKBjzQvQmNEmK9jiYiZUQFXs7sOHCCUbO38EHmDs7m5XP9ZXV5qEcTLq0X53U0EfEzFXg5lXP0FGPnbmP8wiyOnT5LjxaJDOvZlLSUml5HExE/UYGXc4dP5vLegm2MnbeNA8fP0KlxPI9d1YyOjWvqoiCREFdcgV/wKhEzG2tmOWa26pxtNc1supltLPyzhr8Dy8WJi45k2JXNmPsvPfnXX7Vi095j9B29kDvfXMjcjfsI5C9qEQmMklzmNw7ofd623wEznHPNgBmFjyUIxERFMKRbY+b8tif/ftOlbD9wgv5vLeKWEfP5bn2OilykHCnREIqZpQCfO+daFz5eD/Rwzu0ys7rATOdciwt9Hw2hBN7ps3lMztzJiJmbyT50kssbxPHYlc24qlUtDa2IhAifxsCLKPBDzrnqhZ8bcPCnx0W8digwFCApKSk1KyurlP8J4oszZ/OZ8v1OXvtuEzsOnOTSetV49MpmXHtJba0UJBLkyqzACx8fdM5dcBxcZ+Dey83L59PlP/L6d5vYuu84LetUZdiVTbmudV3CVeQiQanUb2IWY0/h0AmFf+b4Ek4CJzI8jNtSGzB9eAYv3dmW3Lx8hr3/Pb1ems2ny7PJy9cYuUioKG2BTwUGFn4+EPjUP3EkUCLCw+jTrj7ThnfntbvbEW7G45OW0+ul2Xy+8kfyVeQiQe+CQyhmNhHoASQAe4A/Ap8AHwJJQBZwh3PuwIV2piGU4JWf7/h69W5enL6BjTnHaFG7KsOvaUavS+vozU4Rj+lCHimRvHzH5yt/5OVvN7Jl33EurVeN4Vc316wVEQ/5ewxcyqnwMOPmtvWZNjyDv9/RhmOnzzLk3Uz6vD6PmZpHLhJUdAYuvyg3L58py7J5ecZGsg+dJDW5Bk9e05zOWrdTJGA0hCI+OXM2n8lLd/DaPzex6/Ap0hvV5MlrmtOxcbzX0UTKPRW4+MXps3lMWryD17/bRM7R03RpGs+T1zQnNVl3PxQpKypw8atTuXlMWLSdETM3se/YGa5qWYvf9GpBq7rVvI4mUu6owKVMnDhzlrfnbePNWZs5evosN7Wpx5PXNCc5XisEifiLClzK1OETuYycvZm3523lbJ7jzvYNeeyqZtSupjU7RXylApeAyDlyilf/uYmJi7cTEW4M7JzCQ92bUD0myutoIiFLBS4BtX3/CV78dgOfLM+mSlQED3RvzH1dGhFbKcLraCIhRwUunli3+wh/+2YD367dQ0KVKIb1bErfDklUigj3OppIyFCBi6eWZh3k+W/WsXDLAepXj+aJq5txyxUNdAtbkRLQpfTiqdTkGky8vyPvDkqnZmwUT3+0kutfnqNl3kR8oAKXgDEzMponMnVYF167ux2nzuZx39tL6DdmET/sPOx1PJGQowKXgDMzbri8HtOHd+fZGy9h3e6j3PjaXB6f9D07DpzwOp5IyNAYuHjuyKlc3py1mTFztuIcDOiUzLCeTakRq6mHIqA3MSUE7Dp8khenb+CjpTuJrRTBIz2bcm/nFCpHasaKVGx6E1OCXt24aP77tjZ89XgG7VNq8l9frePKv83k46U7tVanSBFU4BJ0WtSpyth72/P+/R2Ir1KJpyav4IZX5zJ7w16vo4kEFRW4BK3OTRL49JEuvNK3HcdO53LP2MUMHLuYjXuOeh1NJCiowCWohYUZN7Wpx7dPdudff9WKZdsP0vvlOfzbJ6s4cPyM1/FEPKUCl5BQKSKcId0aM+vpnvTrkMT7i7fT/fnvGD17C6fP5nkdT8QTKnAJKTVjo/jTza35+vFupCbX4D++XMu1L87m61W7dUWnVDgqcAlJzWpXZdx96bwzKJ2o8DAeHL+Uu0YtZFW2ruiUikMFLiGte/NEvnq8G3/u05qNOce48bW5PD15BTlHTnkdTaTM+VTgZva4ma0ys9Vm9oSfMolclIjwMAZ0TOa73/Tg/m6N+WR5Nj3+NpNXZ2zk5BmNj0v5VeoCN7PWwP1AOtAGuMHMmvormMjFiouO5JnrW/Htk93JaJbIC9M3cPXfZ/HlD7s0Pi7lki9n4K2ARc65E865s8As4Bb/xBIpveT4WEYOSGXS0I5UrRzBwxOW0W/MItbv1vxxKV98KfBVQDczizezGOB6oOH5TzKzoWaWaWaZe/fqSjoJnI6N4/n80a78uU9r1uw6wvWvzOHZqas5fCLX62gifuHTzazMbDDwMHAcWA2cds49UdzzdTMr8crB42d4Yfp63l+0neoxUTzdqwV3pDXUikASEsrkZlbOubecc6nOuQzgILDBl+8nUlZqxEbxXJ/L+OzRrjRNrML/+8cP3Pz6XJZmHfA6mkip+ToLpVbhn0kUjH+/749QImXl0npxfPBAR17p2459R89w64gFDP9gOXs07VBCUISPr//YzOKBXOAR59wh3yOJlC2zgvurXNWyFm/M3MTo2VuZtno3j17VjPu6pFApQvcfl9CgBR2kwsvaf5w/f76Wb9fuoXFCLP9+86V0a5bodSyR/6EFHUSKkRwfy5iBaYy7rz35zjHgrcUMe38Zuw9rWEWCmwpcpFCPFrX4+okMhl/dnGlr9nDVCzMZM2cLZ/PyvY4mUiQVuMg5KkeG8/jVzZg+PIP2jWry3BdrueHVuWRu02wVCT4qcJEiJMfH8va97RnZP5UjJ3O5beQCnp68gv3HTnsdTeR/qMBFimFm9G5dh2+f6s6D3Zsw5ftsrnxhFu8v2k6+FlmWIKACF7mAmKgIfnddS756vBst61TlmSk/8OsR83XvcfGcClykhJrVrsqkoR158c42ZB88wU2vzeXZqas5dvqs19GkglKBi1wEM+PX7Row46ke9OuQzDsLtnHN32cxbfVur6NJBaQCFymFuOhI/tynNR8/1JlqlSMZ+t5SHngvU3PHJaBU4CI+uCKpBp8/1pXf9m7BzPV7ufrvs3h3wTby9CanBIAKXMRHkeFhPNyjKdOGZ9C2YXX+8Olqbhs5n3W7j3gdTco5FbiInyTHx/Le4HRevLMNWftPcMMrc/nvr9dxKlfrckrZUIGL+NFPb3J++2R3+rSrzxszN9PrpdnM27TP62hSDqnARcpAzdgo/nZ7G94f0gED+o1ZxJMfLOfg8TNeR5NyRAUuUoY6N03g6ycyGNazKVNX/Mg1L87iyx92eR1LygkVuEgZqxwZzm96tWDqsK7UiavMwxOW8dD4peQc1ZRD8Y0KXCRALqlXjU8e7sJve7dgxrocrvn7bD5eupNALqoi5YsKXCSAIgqnHH75WDea1qrCU5NXcN+4Jfx46KTX0SQEqcBFPNC0VhU+fKATf7zxEhZtOcC1L85mwqIs3eVQLooKXMQj4WHGfV0a8c0TGVzeII7fT1nF3WMWkrX/uNfRJESowEU8lhQfw4QhHfjPWy5jdfYRer00mzFztuhyfLkgFbhIEDAz+qYnMe3JDDo3SeC5L9bSd5TOxuWXqcBFgkjduGjeGpjG325vw9pdR7ju5Tm8tzBLM1WkSD4VuJkNN7PVZrbKzCaaWWV/BROpqMyM21Ib8M3wDFKTa/Bvn6zinrGLNVNFfqbUBW5m9YHHgDTnXGsgHLjLX8FEKrp61aN5d1A6z/VpzdKsg/R6cTaTM3fobFz+h69DKBFAtJlFADHAj75HEpGfmBn9Oybz9eMZtKpbjac/Wsn972bqKk4BfChw51w28DdgO7ALOOycm3b+88xsqJllmlnm3r17S59UpAJLio9h0tCO/OuvWjF74z6ufXE2n63Q+VJF58sQSg3gZqARUA+INbP+5z/POTfKOZfmnEtLTEwsfVKRCi4szBjSrTFfPtaN5PhYHp34PY9MWMYB3eGwwvJlCOVqYKtzbq9zLhf4B9DZP7FEpDhNa1Xh4wc78XSvFkxbs5teL81m1gb967Yi8qXAtwMdzSzGzAy4Cljrn1gi8ksiwsN4pGdTPnmkC9WjIxk4djHPTl2t1X8qGF/GwBcBHwHLgB8Kv9coP+USkRK4tF4cnz3alfu6pDBu/jZufHUuq3887HUsCRAL5JSktLQ0l5mZGbD9iVQkszbs5TeTV3DoxBme7tWCIV0bExZmXscSPzCzpc65tPO360pMkXKie/NEvnkigytb1uIvX66j35hFuvinnFOBi5QjNWOjGNk/lf++9XJW7DxE75dmM1XTDcstFbhIOWNm3NG+IV893o0mtarw2MTvGf7Bco6cyvU6mviZClyknEqOj2XyA5144upmTF3xI9e9NIelWQe9jiV+pAIXKcciwsN44urmTH6wE2Zwx5sLeP27TVr5p5xQgYtUAFck1eCLx7rRu3Udnv9mPQPGLiLniO6nEupU4CIVRFx0JK/1bcd/3XIZS7MOct3Lc5i5PsfrWOIDFbhIBWJm3JWexGfDupJYtRL3vr2E//hiDWfO5nsdTUpBBS5SATWrXZVPHunCgI7JjJ6zldtGzmfbPi3fFmpU4CIVVOXIcP7cpzUj+6eybd9xbnh1Lp8uz/Y6llwEFbhIBde7dR2+eiKDlnWq8vik5Tw9eQUnz+imWKFABS4i1K8ezaShHXn0yqZ8tGwnv35jHpv3HvM6llyAClxEgII5409d24K3723PniOnuOnVuVr1J8ipwEXk/+jRohZfPNaNlnWr8ejE7/nDp6s4fVZDKsFIBS4iP1OvcEjl/m6NeHdBFrePXMCOAye8jiXnUYGLSJEiw8P4/a8uYdSAVLbuO86vXpnD9DV7vI4l51CBi8gvuvbSOnzxaMFCyve/m8l/frmW3Dxd+BMMVOAickFJ8TFMfrATAzom8+bsLdw9eiF7dC8Vz6nARaREfrrw5+W72rL6xyPc8Opclmw74HWsCk0FLiIX5ea29fnkkS5UqRRB31ELGTdvK4FcW1f+lwpcRC5a89pV+XRYF3q0qMWzn63hqQ919aYXVOAiUirVKkcyakAqT13TnCnLs7l1xHxNNQwwFbiIlFpYmPHoVc0YO7A9Ow+e4IZX5zJrw16vY1UYpS5wM2thZsvP+ThiZk/4MZuIhIieLWvx2aNdqRtXmXvfXqxl2wKk1AXunFvvnGvrnGsLpAIngCn+CiYioSU5PpZ/PNyZm9rU4/lv1vPg+KUcPZXrdaxyzV9DKFcBm51zWX76fiISgmKiInjpzrb84YZLmLEuh1vemE/Wfi0UUVb8VeB3ARP99L1EJISZGYO6NuK9QensPXaam16bx7xN+7yOVS75XOBmFgXcBEwu5utDzSzTzDL37tWbGyIVReemCXz6SBdqVa3EPWMX8878bZov7mf+OAO/DljmnCvyLjfOuVHOuTTnXFpiYqIfdicioeKncfGeLRL549TVPDNllRZQ9iN/FHhfNHwiIsWoWjmSUQPSeLhHEyYu3k7/MYvYf+y017HKBZ8K3MxigWuAf/gnjoiUR2Fhxm97t+Tlu9qyYuchbnptHmt3HfE6VsjzqcCdc8edc/HOucP+CiQi5dfNbesz+cFO5OU7bh0xn69X7fI6UkjTlZgiElCXN6jO1GFdaF67Kg+OX8br323Sm5ulpAIXkYCrVa0yk4Z25Oa2BRf9/MvHK/XmZilEeB1ARCqmypHhvHRnW5LjY3llxkZ2HjzJiH6pxMVEeh0tZOgMXEQ8Y2Y8eU1zXri9DUu2HeCWEfPYvl93NCwpFbiIeO7W1AaMH9yB/cfP0OeNeSzN0ko/JaECF5Gg0KFxPFMe7kK1yhH0Hb2IT5dnex0p6KnARSRoNEqIZcrDXWjboDqPT1rOqzM2aobKL1CBi0hQqREbxXtD0vl1u/q8MH0DT3+0ktw8zVApimahiEjQqRQRzt/vaENyfAwvfbuRnKOnGdHvCmIrqbLOpTNwEQlKZsYTVzfnr7dexrxN+7hz1AJyjp7yOlZQUYGLSFC7s30So+9JZXPOcW4dMZ8te495HSloqMBFJOhd2bI2E4d25PjpPG4dMZ9l2w96HSkoqMBFJCS0bVidfzzUmWrRkdw9eiHfrilyCYIKRQUuIiEjJSGWjx/qTPPaVRn6XibvL9rudSRPqcBFJKQkVKnExPs7ktE8kWem/MCL0zdU2LniKnARCTmxlSIYfU8at6U24OUZG3l26mry8yteiWtSpYiEpMjwMJ6/7XJqxEQyes5WDp/M5fnb2xAZXnHOS1XgIhKyzIxnrm9F9Zgonv9mPUdOneWNfldQOTLc62gBUXF+VYlIuWRmPNKzKc/1ac1363O4Z+xijpzK9TpWQKjARaRc6N8xmVfuaseyrIP0HbWQfRVg5XsVuIiUGze2qceYgWls3nuMO0YuYOfB8r04hApcRMqVHi1qMX5wB/YdO83tIxewKaf8XnqvAheRcictpSYfPNCJ3DzHXaMWsG73Ea8jlQkVuIiUS63qVuPDBzoSERbGXaMWsir7sNeR/M6nAjez6mb2kZmtM7O1ZtbJX8FERHzVOLEKHz7QidioCPqOXsj35ewmWL6egb8MfO2cawm0Adb6HklExH+S4mP48MFO1IyNov+YRSzeWn4WTC51gZtZHJABvAXgnDvjnDvkp1wiIn5Tv3o0Hz7QiTpxlRk4djHzNu3zOpJf+HIG3gjYC7xtZt+b2Rgziz3/SWY21MwyzSxz7969PuxORKT0alerzKShnUiOj+G+cUv4bn2O15F85kuBRwBXACOcc+2A48Dvzn+Sc26Ucy7NOZeWmJjow+5ERHyTWLXgTobNa1dh6LuZTFu92+tIPvGlwHcCO51ziwoff0RBoYuIBK0asVFMGNKRS+vF8fCEZSFd4qUucOfcbmCHmbUo3HQVsMYvqUREylBcdCTvDk6ndf04Hnl/Wciu7uPrLJRHgQlmthJoC/zF50QiIgFQrXJBiV9StxoPTVjKP9eFXon7VODOueWF49uXO+f6OOfK1yRLESnXCkq8Ay3rVOPB95aF3BubuhJTRCq0uOhIxg/uQPM6VXjgvaXM2hA6s+VU4CJS4cXFFJR408Qq3P9uJnM2hkaJq8BFRIDqMVFMGNKBxgmxDHknMyQu9lGBi4gUqhEbxfv3d6RRQiyD31nCkm3Bfdm9ClxE5Bw1Y6MYP6QD9apHM+jtJfywM3jvYqgCFxE5T0KVSowf3IFq0ZHcM3YRG/cc9TpSkVTgIiJFqFc9mglDOhARHka/MYvI2n/c60g/owIXESlGSkIsE4Z04ExePv3GLGLX4ZNeR/o/VOAiIr+gee2qvDsonUMncuk/ZlFQrXavAhcRuYDLG1Rn7L3tyT50knveWszhk7leRwJU4CIiJZLeqCZvDkhjY85RBo1bwskzeV5HUoGLiJRU9+aJvHxXO5ZtP8ijE5dxNi/f0zwqcBGRi3D9ZXX5002X8u3aHH4/ZRXOOc+yRHi2ZxGREDWgUwo5R0/z6j83UataJZ66tsWFX1QGVOAiIqXw5DXN2VtY4glVKjGwc0rAM6jARURKwcx4rk9r9h07w7OfrSahSiV+dXndgGbQGLiISClFhIfx2t3tSEuuwfAPljN/c2DvYKgCFxHxQeXIcMbc056UhBiGvruUtbuOBGzfKnARER/FxUTyzqB0qlSKYNC4Jew5ciog+1WBi4j4Qd24aMbe254jJ3MZNG4Jx0+fLfN9qsBFRPzkknrVeO3uK1i76wiPTfyevPyynSOuAhcR8aOeLWvx7ze3Zsa6HP702eoyvdBH0whFRPxsQMdktu8/zug5W0mKj2Vw10Zlsh+fCtzMtgFHgTzgrHMuzR+hRERC3f+7rhU7DpzkuS/W0LBGNNdeWsfv+/DHEEpP51xblbeIyP8KCzNevLMtlzeozmOTvmfFjkP+34ffv6OIiAAQHRXOmHvSaJ9Sk6qV/T9ibb4MsJvZVuAg4IA3nXOjinjOUGAoQFJSUmpWVlap9yciUhGZ2dKiRjl8PQPv6py7ArgOeMTMMs5/gnNulHMuzTmXlpiY6OPuRETkJz4VuHMuu/DPHGAKkO6PUCIicmGlLnAzizWzqj99DlwLrPJXMBER+WW+jKrXBqaY2U/f533n3Nd+SSUiIhdU6gJ3zm0B2vgxi4iIXARNIxQRCVEqcBGREKUCFxEJUT5dyHPROzPbC/h6JU8CENh1iy6eMvqHMvqHMvqHlxmTnXM/u5AmoAXuD2aWGez3XVFG/1BG/1BG/wjGjBpCEREJUSpwEZEQFYoF/rMbZgUhZfQPZfQPZfSPoMsYcmPgIiJSIBTPwEVEBBW4iEjICpoCN7OxZpZjZqvO2Xa7ma02s3wzK3b6jpltM7MfzGy5mWUGOOPzZrbOzFaa2RQzq17Ma3ub2Xoz22RmvwvSjF4exz8X5ltuZtPMrF4xrx1oZhsLPwYGaca8wucsN7Opgcx4zteeMjNnZgnFvNaz43gRGT07jmb2rJlln7P/64t5bUB+rovlnAuKDyADuAJYdc62VkALYCaQ9guv3QYkeJTxWiCi8PO/An8t4nXhwGagMRAFrAAuCaaMQXAcq53z+WPAyCJeVxPYUvhnjcLPawRTxsKvHSvrY1hcxsLtDYFvKLho7mf/P70+jiXJ6PVxBJ4FfnOB1wXs57q4j6A5A3fOzQYOnLdtrXNuvUeRfqaYjNOcc2cLHy4EGhTx0nRgk3Nui3PuDDAJuDnIMgZMMRmPnPMwloJl+s7XC5junDvgnDsITAd6B1nGgCkqY6EXgd9SfD5Pj2MJMwbML2S8kID9XBcnaArcRw6YZmZLC9fg9Mog4KsittcHdpzzeGfhNi8UlxE8Po5m9h9mtgPoB/yhiKd4fhxLkBGgspllmtlCM+sTuHRgZjcD2c65Fb/wNE+PYwkzgofHsdCwwiGzsWZWo4ive/73sbwU+AXX5ixrZvZ74CwwIdD7LqkSZPT0ODrnfu+ca0hBvmGB3HdJlTBjsiu45Ppu4CUzaxKIbGYWAzxD8b9YPHeRGT05joVGAE2AtsAu4IUA7rvEykWBO4/X5jSze4EbgH6ucHDsPNkUjPn9pEHhtoApQUbPj+M5JgC3FrHd8+N4juIynnsct1Dw/k27AGVqAjQCVpjZNgqOzzIzq3Pe87w8jiXN6OVxxDm3xzmX55zLB0ZT9M+C538fQ77AzeO1Oc2sNwVjeTc5504U87QlQDMza2RmUcBdQJm9q16ajEFwHJud8/BmYF0RT/sGuNbMahT+k/bawm0BUZKMhdkqFX6eAHQB1gQin3PuB+dcLedcinMuhYJ/0l/hnNt93lM9O44lzejlcSzcZ91zHv6aon8WPP25BoJqFspECv6pkkvB/9TBFBy4ncBpYA/wTeFz6wFfFn7emIJ3f1cAq4HfBzjjJgrGwZYXfow8P2Ph4+uBDRS8ax10GYPgOH5MwQ/JSuAzoH7hc9OAMee8dlDhf88m4L5gywh0Bn4oPI4/AIMDmfG8r2+jcIZHMB3HkmT0+jgC7xXudyUFpVz3/J+ZwscB+bku7kOX0ouIhKiQH0IREamoVOAiIiFKBS4iEqJU4CIiIUoFLiISolTgIiIhSgUuIhKi/j/XIpX9OX/xNwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "theta = np.arange(-5, -4, 0.0001)\n",
    "r = [r_t(q) for q in theta]\n",
    "x = [x_t(q) for q in theta]\n",
    "y = [y_t(q) for q in theta]\n",
    "\n",
    "fig1 = plt.figure(1)\n",
    "axes1 = fig1.subplots(1, 1)\n",
    "axes1.plot(x, y,label='parametric curve')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ec8b644",
   "metadata": {},
   "source": [
    "Генерация наблюдений "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "73b0f570",
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.arange(-5, -4, T)\n",
    "Y = []\n",
    "for th in theta:\n",
    "    rr = r_t(th) +  np.random.normal(0, sigma_r)\n",
    "    xx = x_t(th) +  np.random.normal(0, sigma_x)\n",
    "    yy = y_t(th) +  np.random.normal(0, sigma_y)\n",
    "    Y.append(np.array([rr,xx,yy]))\n",
    "axes1.plot([j[1] for j in Y], [j[2] for j in Y], 'o', color='black',ms=1);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "dc2b478c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjrElEQVR4nO3deXxU5dUH8N+DuEZcIIgLkqhg3UWJrUuxValitNG6tFqtMC7U2ljq6wthUYdURUlqscqr1oWIu3VDGgIIImARxbCKJSwqq0oyFkTCmsx5/8gMTCaz3vW5c3/fz2c+JJM79565JOeeOfd57lUiAiIi8p52bgdARETGMIETEXkUEzgRkUcxgRMReRQTOBGRR7V3cmP5+flSWFjo5CaJiDxv3rx5IRHpHP+8owm8sLAQtbW1Tm6SiMjzlFKrEz3PFgoRkUcxgRMReRQTOBGRRzGBExF5FBM4EZFHMYETEXkUEzgRkUflTAIPhUKorKxEKBRyOxQiIkfkTAKvqqrC4MGDUVVV5XYoRESOcHQmpp0CgUCrf4mIcl3OJPD8/HwMGjTI7TCIiByTMy0U0g/PSxDZiwmcbMPzEkT2ypkWCumH5yWI7MUKnGwTPS+Rn59vyfrYkiFqjQmc0tIlcbIlQ9QaWyjURigUQlVVFQKBAPLz83cnTgCujvRhS4aoNVbgNtGlajUivtINBAKoqKhwPXFa3ZJxipd/F0hvrMBtokvVakR8pZvNGPv46p28/btAemMCt4kOH/eNJlMzk6KYrNrS4XeBclPaBK6UGgvgcgD1InJK5LlKAL8EsBPAFwACIrLJxjg9R4eZoW4kUyartnT4XaDcpEQk9QJKnQ9gC4AXYhL4xQCmi0iTUmoUAIhIWbqNFRUVCe9K75xoBV5SUoIJEyawrUHkUUqpeSJSFP982pOYIjILwH/jnntPRJoi334MoKslUZKlopXfhAkTOPyOKAdZ0QO/GcDryX6olBoAYAAAdOvWzYLNUbYyaWvw5COR95gaRqiUGg6gCcDLyZYRkadFpEhEijp37mxmc57l9jCyTIbfcZIMkfcYrsCVUv3RcnLzIknXSPc5L4zM4MlHIu8xlMCVUn0BDAbwMxHZam1IuccLydGqkRJsxRA5J20LRSn1KoA5AH6klFqnlLoFwBgAHQBMVUotVEo9ZXOcnubVGYRGsBVD5Jy0FbiIXJ/g6edsiIUs4HYF7IVPG0S5gtdCyTFuV8B++rRB5DZOpc8xrICJ/IMVuM2cHkLICpjIP5jAbeZ2S4OIchdbKDbzQkvD7ROfRGQMK3Cb2dnSsKo9w08JRN7ECtzDrJrh6YVPCUTUFitwD7PqVmdWfEpw+3oviegYE5GVmMAjvPjHrtOIEx3bMDrGRGQltlAivHDBKZ3p1oYJhUJobGxEMBjUJiYiqzGBR+iWgKy2KtSIr0KNuOCEw2xZv263DauqqkJ5eTkqKiq0+IRCZAcm8AjdEpDVnv33l3jp4zX42fGdcc9lJ6JHlw5uh5SWmeGNuX5AJgLYA/eN+y4/GfdcdiLmr9mIvn//EMF3l2Bj405HY4ieZ1i2bFlG5xvM9LB1Oj9AZBdW4D6xT/t2uLX3sbjqzK4YPXU5Xvx4NcYv/Bp/7tMDN55dgL33sv9YHk3IM2bMQE1NDYDU5xtYRROllvau9FbiXen1sezbH3B/9X8wc/EX2PfLWfjr8IG48pwTky5vxWzN6DpKSkowYcIEX8785KxXMiLZXekhIo49evXqJaSPcDgst919rwCQQ34ekN8994msrP8h4bIVFRUCQCoqKhyO0hoNDQ1SUVEhDQ0Nrm7L6/uR3AGgVhLkVLZQfEwphZFD/oRjOh+I/U66EGNrv0PfR2fhlp8eizsv7I68fff8emTaznCzwozfduz3Tg4TTbUttoXIUomyul0PVuB6q9+8Xe7+50IpKKuWnzw4TSYsXC/hcDirdVhRYRqtluO3Hfu9VRV4Jutxstonf0CSCpwJXCO6/OHXrvpOiv8+SwrKquW6f8yRZd9uzvi12byHZMsmOggkWjb+uXTfW4EtEHIDE7gHWJUcrEhcTc1heWHOKjltxBQ5buhEuf9fn8vmbTtNxRUv2ftNFH+iZd1IprocZMlfmMAtYucfsFXrtjKxfbdlhwx5a5EUDqmWogemytvz12bdVkkm1fvNpJpmMiW/YAK3iB1Vn9WJyI7EtmDNRil5/EMpKKuWXz/1kazYkHi0SqbSxchWBdEeTOAW0b2vamdV2twclpc/Xi2nBidL92ET5ZH3lsm2nU2G1pXuPbO6dh73ub6YwB1m9GSe2T8ioweDbGKo37xdBr46XwrKquXnlR/Iv1dkH7fOyULn2OzETz36SpbAOQ7cJtmMO469kFZlZaWp8cpGxxnHxgsgZQydO+yLR687A1f36op7xi/BDc9+gm7rpuHDlx81HLdO/HppYY5R96BEWd2uBytw+15nVux26+rqpLi4WOrq6tK+btvOJvnrlDo55s+vyBG/uFWefm+BNDenP8lpZ7Vn9hONXytw0hfYQrGGH/64jSTX5d9ulmuenC0FZdVy7ZMfyYoNyceONzQ0SDAYlGAwaMt+jI3fD20BP/xO+h0TuEWYEJJrbg7La3NXy2kjpkiPYTUyZvoK2dnU3Ga5bPah2Qpah+Rmdwx++J30OyZwi+iQEKxg50nH+s3b5Y6X5klBWbVc+ugsWbJ+k+F1GUlOuv0f2Z1gdXu/ZD0mcGol26RiJAlN+uxr6XX/VDlu6ESpnFwn23dlNuTQbAVtddI3myCZYMksJnAylBijy9XV1RlKQhsbd8hdry+QgrJqueiRGTJv9X9tn8RjddJni4Lc5ukE7pcKRsdeqdFqNv4k5fS6DXLOyGlSOKRaLg7c7egknkzWZ2cFbjY2IsMJHMBYAPUAlsQ81xHAVAArIv8emm49YiKB+6UC0rFXaqaajX8vm7ftlOHvLJaud74sxxQPkCm1y01tM/ZAkeoTQjSe6LJ2jX4xwi+/22SOmQR+PoAz4xJ4BYAhka+HABiVbj3CCjwto7M33ZIshnTDBD9aGZLeo6ZL4ZBqeaD68zbT8TNNarEHij59+uxO0sniDAaDCQ8sZit0M3T4fyT9mWqhACiMS+DLABwR+foIAMsyWQ974MYk+iPXoXIzE8OW7btk2NuLpaCsWvo8MkMWr920+2dGKvBBgwYlTeCJls92X8ZW8TokXCZ+f7E6gW+K+VrFfp/gtQMA1AKo7datm3PvOIckSjA6/AFbEcOMZfXy4wdbRqo8OnV5wnHjdseSTQUereLdbnnocAAn59iWwCPfb8xkPazAjdEhWZuV6j1saty5++JYv3z8w6SzOHXYDzrEoFMc5Ay2UMhVmVSMExd/LT3Lp8jxw2vk2Q+/bHNNFbuvn6LbCU4rMNHnhmQJvF38xa0yNAFAv8jX/QC8a3A9OSsUCqGyshKhUMjtULQQCARQUVGBQCCQcN+EQiF8PvklvHLTSejdIx/3V/8Hv332Y3y9aVvCdSSSzT6PX7aqqgrl5eUoLy9HVVWVyXerj+iVFXPpPVGMRFldWlffrwL4BsAuAOsA3AKgE4D30TKMcBqAjunWIz6rwNmjTC7Rvol9LhwOy+ufrpET750kp42YIjWLvza8XpHMTgKzAiedwcsTebyIfzjJJdo3iZ77qmHL7tu4lb25SBp37Mp6vSL6ngQmylSyBK5afuaMoqIiqa2tdWx75H27msMYPXU5npz5BY7plIfHrj8Dpxx1cMrXhEIhVFVVIRAIID8/v833dnBiG+RfSql5IlIU/7zRHjiRI/beqx0G9z0Br9x6NrbubMavnpiNp2d9gXA4eeER3/eN3vHIzsTKXjO5gQmcPOGc4zph8p9746ITumBkTR1uGjsXGzZvT7hsupOdmcr0pGgoFEJjYyOCwSBvR0aOYgL3Ea+PjDnkgH3w5I1n4uGrTsW81RvR99FZ+GBZfZvlrKq4M62qoyNY8vLy2D4hRzGB+4hbH/OTHTiMHFCUUrjux93wrzt/ii4H7YdA1aeomFyHpuaw4TiSybSSt6riJ8paojObdj38NArFLWYui2rXyIxkw/vMDrXctrNJhrzVcj2Va56cLd9s2mYoDiLdgcMI/cFMksrmtVZcOdGqA8Y789fJifdOkjP+8p58ULfBVMxGDnJuDUnkUEj/YAL3MKcuM5vNa3WrZlfW/yCXjJ4pBWXVMmrSUtll8KJY6d5XuklIdtH1ipTkDCZwGzhVAUWvgJfqUqlOc6sdk0pLS2VR2paK1W0mJ94rJyP5GxO4DZyqgHRM4Om4WR1GWypn/uU9mb0y+R16Mo3NyoOV0aTLZO1vTOA2cOqPyot/vG7G3NDQIIPvu19+Wj5ejhlSLf+YuVLC4XCrn2cTWzThFxcXJ3xNNgcEqw9sXvzdoOwxgZM27E460SR5/8iH5fYXa6WgrFrueGmebNme+loqyTQ0NEhxcXHSxOvUOYpE2Af3ByZw0oaZpJPtSJJwOCyV4+fKoRcE5Pzy8fJF/Q+GtmNl4rVyNEzszZxZjecuJnDShplEYyT5R19z5MW3yin3TZYpS77JejtWJsdM3oOuo2HIHckSeHsb5wiRpty+cl50qnsqya4oWFJSgsbGRjQ2NiIUCmUUf3SGZN+rrsPwSasx4MV5KL2gO+76xfHYq51K+Zrov9FZrADSxp5pPKlmbqZbJtHPM1kv5ZhEWd2uBytwPXihUouPMfb72K8zrYyjy639+lsZ/EbLUMP+Yz+Rzdt2ZhQP2xPkJrCFQlFeSEapetCxX2d6MIpf7sU5q+S4oRPlokdmyFcNWxx5D27TLR7KHBM4mWLXH7/Z9WZbgccuN3tlg5xePkVOL58is1dYn9SMfNKxej8bOdiRfpjAfc5sYrBr/HJ0klKq9do563NVaIv0eWSGHDt0orzw0Vem15fsk0KmrN7PRtpNuSDX3isTuM+Z/UO2a/xyMBhMu14jIzKysXnbTrm5aq4UlFXLsLcXy0MPj3LkgmCJpJqqHztk0Mz6/CDXPm0wgfucHR+lnbpwVkND6jvGW5GkmprD8lDNUikoq5YrH6mREQ+MNPy+rL67fexM0FxKSnbKtQMXEzjtZtUvt5NVjt3biu6T599fJD2G18j5FdMznvQTz652k5EKnHJDsgTOu9KTYU6OJ7d7W5WVlRg8eDAqKipw4a9vxW0v1CIsgmduKsJZhR2zjnXMmDEAgNLS0lbxuj0Gn7wp2V3pWYFTTrB6NMuq0Ba5oPID6TGsRsYvWJf1OpJV4bnWmyVngDMxKZeZnSkZPzu0oFMe3r7jXPR74n30v+teDBl4OwZfcRaUSjxzMz6GQCCQcMYoZ0uSlXhT4xzh1B3nnbyzfTbbSnZj4VTrSLf+Qw7YB2c3L8amGVV4+O9PYchbn2FXipsnx8aQn5+PvLw8lJeXt7qJdPRAwfYJWSJRWW7Xgy0U+zj10VyXE5eZtkxSrSPZz+LHco8aNUrK/zlHCsqq5cZnP5bvLZ5+n2sjJsh64CiU3OZUEnAy2aTaVqYHklQjOJKtP9m6X/90jRw3dKJcMnqmfP7l2qRT/bNlxYHKCTrF4jdM4JRTsk0m2XxySLXuWcvr5aR7J0nBpQNarS86o3TQoEFJX5tsvZkeqNxOoDwB6x4mcPI1K5PforUb5bSyN+Soi2+VGYtWisieBN6nT5+kSS52mUzj0OlaJm4fQPwsWQLnKBTyhdhRJmbHYp/W9RC8c3df3DS2I+54cwWe2v8glJaWIi8vDyUlJZgwYULKUSbTpk3bvf10ccTGzREs1EairG7XgxU46SCb67Ck8u332+SS0TOl+7CJ8u7C9WmXj59m73ZFnS2vxZtLYEcFrpS6C8CtAATAZwACIrLd7EGFyE7RCraxsdHU2PEuB+2H139/Dm4bV4uBry3Af7fsQP/zjkm6fH5+PkaMGNEmDq9U1Kni5QxTlyTK6pk8ABwF4CsA+0e+/yeA/qlewwqcdGJVT3fbzia5bdynUlBWLY+8t0zC4bBW8TmB1bm9kKQCNzuRpz2A/ZVS7QEcAOBrk+ujHObkJKBMxPaX08UVH3vs9/vtvReeuOFMXNurKx57fwVG1iyNFjmmRGd2xk4EypTT+zrZRCqyWaKsnukDwEAAWwA0AHg5yTIDANQCqO3WrZtTByyyiJVVoK5VWiZxxS+T6DXNzWG5b/xnu68r3txsrhK3a2w5eQ+sHkYI4FAA0wF0BrA3gPEAbkz1GrZQvMfKRKBrSyCTuOKXSfaacDgsD09qua74Xa8tkF1NzZbH4uR6yBpm/z/sSODXAngu5vubADyR6jVM4N6jUyLI9iYQbsb9+PvLpeudL8vZ1/1J1n+zIePY/FQ5u/1/5CSz/692JPCfAPgcLb1vBWAcgDtTvYYJnMzI5o9Ah0R4zR+GCAApurZUtu1syig2JrXcpF0F3rJOlAOoA7AEwIsA9k21PBN47tHl2ihuxpUqhutLh8qRtz0lPa/+o6xZ/62jsemwD1LRPT6d2JLAs30wgeceP1VRsVLdfDg+Id00cLgAkNOv+qNs3dEUvyrb+PX/JhclS+C8HjiZosPwMTeGJyYa4hf/XDSuYX/4HW4aOBwbu56LW8Z9im07m22LK3ZfWPF/k82+1W2YqC8kyup2PViBkx3cqDQzqcDj43p7/lopHFItv31mjm2VuNX7wmvnHXIV2EKhXKVrLzX+xhAVFRVSNW2RFA6plhue+bjViU07tmnF8l4775CrmMCJXBRbnb5Zu7ZVEncz8UUvcRsMBh3fNmUuWQJnD9yn/NCv1Ok9xvajr+7VFZXXnI7ZX4Rw2wu1ePrZ5wxPmXeCTvuRWmMC9ykz19nQVXyi0ek9xt/M+JpeXVFx9Wn4cEUISzv0wj333rf7DvaZsiKxlpaWoqKiAqWlpUmX0Wk/UpxEZbldD7ZQ9JGL/cr4k2heeI/jPvpKCsqq5bzfDhQAUlxcbMtt4szwwn7MdUjSQlEtP3NGUVGR1NbWOrY98hcvXpM6FArhtnsq8UnzsThg/sv4Yv6HCAaDyMvLS/s+vPh+yRil1DwRKYp/ni0UyhnxbQovqKqqwvh/VOAMWYkd592OvjffDRHJqGWRyftl/zq38Z6YRC6KTrLp378/Hv9oA16YczB6n56PiooDLZkcFe1fA8buOmQ1fmqwFitwogzYVclGq+jOnTtjxC9PxtVndsXTc0PodM41liQ4HWbKxnL6hGiufwJhBU6UgWjiaWxszKg/bUS7dgqjrj4VW3c24S/V/0GnA/fBFT2PMrXO2LsO6SAQCKCxsXH3iBu7q3DdPoFYjQmcEuJH3dbib4RsVyJvv1c7jP5NT2zcOhf/+8YidMrbFz/tkTv7Pz8/H3l5eRg8eDDy8vJsT6peu3F01hINTbHrwWGE3sHrWiQWHVIXncGYbv8YnYr+/badcsnomXLSvZPks3WbrArfEXV1dVJcXCx1dXUJf85hidkDp9JTNvhHllqm+yd6IAwGg2mXjz9ofrNpm5z70Pty+pA3ZWjwgaxu+eam4uLi3WPayRpM4ORpOiWobGRTsSd6jys2bJYjfnFry0Hg/pFJX6vTJ6Z0FThljwmcPKuhoWF3VRefoLyS2M3EOW3+Cul80c1yyUP/ki3bd1m+ftJfsgTOYYSkvaqqKtTU1KC4uLjNySivXKfDzCSji87ojlceG4nl3ysMfG0BmsPOzZ4mvTGBk/aiY5nHjRvXJgHqNs7ZLn1O6oLykpMxbWk9RtYsbfNzJw5kdoypNrvOXB/nnVaistyuB1so5JZcaTGMmLBECsqq5YU5q1o978T7s6PPbnadOvX+7YQkLRSOAydfyJUJHfdcdhJWf7cVIyZ8jm4dD8DPju8MwJkJO3aMqTa7zpwf550Gr0ZIvpBLE5O27GjCNU9+hPUbt+HNP5yLHx3eIeXyufTe/YpXIyRfc/NKhVb3eQ/ctz3G9j8L++2zF25+/lM0/LAj5eut6I/7vtesKSZwIhtEE96yZcvQr18/Uwk0UQI+8pD98Vy/InzXuAO/f7EWO5qak77eihO9Xhnt4zuJGuN2PXgSk5zi9knL6Mm12FmJRmNJ9V5emfmZHPLzgNw5doaEw2HL159uGbf3s1+AE3nIT9wenRBNbHV1dbYmuOj7POjc6+W6Pw5NuJ10SdbMvsrmUgFkXLIEzlEolJPcHp0QOyrEztEhgUAAYRG8OnsFXvu/h5DfYV88/lCw1TKpRuCEQiE0NjYiGAwa2lfxV2lMtA0n+e6EbaKsbteDFTiRPb5Y87Uce9nv5fQhb8jXm7a2+lmqCjxR9W2kLaJLK8XtT152AafSEzkj2xEbVozwOPboIzD5uUrs2vtA/P7Fedi+a89JzVQjcAKBAILB4O4bLADGTljqcj9Sv8zMjWICJ7JYtgnQSMJMlPR7dOmAv/2mJxav+x73jF/ScpIrjegNFsrLy3dv38tJUJcDiVPYAyfL+a4PGSeb/nu2Pejovm1sbER5eTmA1j3nS04+HH+6sDsem74SPy7siF+fdXTW8Tp9GzYrf19897uXqK9i14M9cH/I1T5kOkb6wNnuq0xGfTQ1h+W3z8yR44fXyOfrvzcVnxOs/H3J1d892DGMEMAhAN4EUAdgKYBzUi3PBO4PuiYKuxlJHtnuq9jlUw1VrN+8Xc56YKr8vPID2bxtp+H4jMbp1rpz9XfPrgQ+DsCtka/3AXBIquWZwMkJbv0RO73d+MlC8Yn54y9CcuzQiXLHS/MkHA6bOljkamXrFZYncAAHA/gKkQtiZfJgAicn+CHZNDQ0SDAYlGAwmHKy0BMfrJSCsmp5fvZXWW8jdj+aOTgZvbEz7WFHAu8JYC6A5wEsAPAsgLwEyw0AUAugtlu3bg6+ZfIrPySBTMdvNzeH5eaqudJ92ERZsGZjVtuwaj9mc0D1w8HXCDsSeBGAJgA/iXz/dwD3p3oNK3Cygx8SdrxE7zlZ8lu5er0UXDpAzr7vHfkhyT017Ywxm8sJ+PH/MhN2JPDDAayK+b43gImpXsMETnbQvWpzKikl2k7sDaEPvSAg//P6QltjiKX7/4uX2HUS80MAP4p8PQJAZarlmcDJDrpXbVYkMqPvMfZEZ/kbc6SgrFreXbjecBzZxGT257SHXQm8Z6S/vRjAeACHplqeCZz8yIpEZfSEYuyyu5qa5Vf/9285JThZ1nzXaCgOK0emsELPnC0JPNsHEziRMVYlzjXfNcrJ902Wq5+YLbuamrN+fWxFb/ZSuazAM8cETpQjzCa+d+avk4Kyanl06nJDY8OTjTu3k9+TfbIEzotZEWkkkysTxl6wKdHy6dZx5RlH4cqeR+Kx6Svw4OgnsrqQVn5+PsaNG5f2YldW30Mz3QW/fHvPzkRZ3a4HK3BygperNaPXRoldPpN1bNq6U37y4DQ5v3y8PPjQw5bvK6v723beVcgLwDvykF+kugON7rK9k1Ci5TNZx8H7742Hrz4V/as+hTqrxPIr91l9R6R0V0h0+w5MblEtyd0ZRUVFUltb69j2yJ+8dElRt2Md+vZivP7pWrxx+7noVXCo49unzCil5olIUfzz7IFTzvHSRf2N3MzBSsOKT8QRB++PQW8sanUXn2Ss6DX7tl9tAyZwIhelu/uN3cmuw357o+Ka0/BlqBGVU5alXd6KA47bB61cwh44kYvS9XZj+/mBQMCWdst53fPxu7MLMHb2V7j0lMNRVNgx6bJmes3RdlFJSYnhdVBrTOBEGotNmHaenB1y6QmYXlePoW9/hol/6o192if+cG7mdmtePrmsK57EJPIIu094Tq/bgH5PvI9eOxbh2QetP4fg9glbL+NJTCKPie9/231y9sITuuDo0Fy8+3QF/vr4U5av30snl72CCZxIU/En+5wYvfHcyP9Flz63YPnBRUj26ZyjSPTBBE6kqfgRKkZHb2STcE8s7Iq/PXAf5tc348156xIuw1Ek+uBJTCJNxZ8wNDoCJNuTh9eddTTeWbAOD9YsxYUnHIZOB+7b6ud+nfWoI1bgRB6Rn5+/ezRKNu2LdGPN47Vrp/DQVafi8IP2Q/0POxLGwV62HpjAiTzESPsitpLPtJXS/bAOmDSwN0484qCUy7Ef7i62UIg8xEz7IlkrJdnwPqWU4XWSM5jAiTzEzESaZMnfTBJmP9xdnMhD5BJdJrboEgclx4k8RBaxqu/r9nC86PsAwJOSHsUWClGWrOr72tV+yLSiZv/a+5jAibJkVeI1089OJdPEzP6197EHTqQxI/1p9rRzD3vgRB5kZtw3k3fuYwuFSGNsc1AqTOBEGrOrT065gS0UohzAKe3+xAROlAPcHlNO7mALhcgmTo4GYa/cn1iBE9lkzJgxGDx4MMaMGWP7tjjyxJ+YwInIU9jv34MtFCKblJaWIi8vj20Ni/ESAHuYTuBKqb0A1AJYLyKXmw+JKDdkMgSQsyazx37/Hla0UAYCWGrBeoh8x++jR4y0Q9jv38NUBa6U6grgMgAPAvgfSyIi8hG/V5Nsh5hjtoXyKIDBADqYD4XIf/w+09LvBzCzDLdQlFKXA6gXkXlplhuglKpVStU2NDQY3RwR5SC2Q8wx0wM/D0CJUmoVgNcAXKiUeil+IRF5WkSKRKSoc+fOJjZHRESxDCdwERkqIl1FpBDAdQCmi8iNlkVGREQpcSIPEZFHWTKRR0RmAJhhxbqIiCgzrMCJiDyKCZyIyKOYwImIPIoJnIjIo5jAiYg8igmciMijmMCJiDyKCZyIyKOYwImIPIoJnIjIo5jAiYg8igmciMijmMCJqA0j96ok5zGBE1Ebfr/ZsldYcjlZIsotvFelNzCBE1Ebfr/ZslewhUJE5FFM4EREHsUETkTkUUzgREQexQRORORRTOBERB7FBE5E5FFKRJzbmFINAFYDyAeg4xxdxpUdxpUdHePSMSaAccUrEJHO8U86msB3b1SpWhEpcnzDaTCu7DCu7OgYl44xAYwrU2yhEBF5FBM4EZFHuZXAn3Zpu+kwruwwruzoGJeOMQGMKyOu9MCJiMg8tlCIiDyKCZyIyKMcT+BKqYFKqSVKqc+VUn92evsxcYxVStUrpZbEPNdRKTVVKbUi8u+hmsR1bWR/hZVSrgxhShJXpVKqTim1WCn1jlLqEE3iuj8S00Kl1HtKqSPdjinmZ3crpUQple9kTMniUkqNUEqtj+yrhUqpYh3iijx/Z+T363OlVIUOcSmlXo/ZV6uUUgudjiuWowlcKXUKgNsA/BjA6QAuV0p1dzKGGM8D6Bv33BAA74tIDwDvR7532vNoG9cSAFcBmOV4NHs8j7ZxTQVwioicBmA5gKFOB4XEcVWKyGki0hNANYD7NIgJSqmjAVwMYI3D8UQ9jwRxARgtIj0jjxqHYwISxKWUugDAFQBOF5GTAfxVh7hE5DfRfQXgLQBvuxDXbk5X4CcC+EREtopIE4CZaElMjhORWQD+G/f0FQDGRb4eB+BKJ2MCEsclIktFZJnTscTFkCiu9yL/jwDwMYCumsS1OebbPACOnqlP8rsFAKMBDHY6nqgUcbkqSVx/APCwiOyILFOvSVwAAKWUAvBrAK86GlQcpxP4EgC9lVKdlFIHACgGcLTDMaTSRUS+iXz9LYAubgbjMTcDmOR2EFFKqQeVUmsB3ADnK/BE8VwBYL2ILHI7lgRKIy2nsW60DZM4Hi254hOl1Eyl1FluBxSnN4ANIrLCzSAcTeAishTAKADvAZgMYCGAZidjyJS0jK/kGMsMKKWGA2gC8LLbsUSJyHARORotMZW6GUukWBkGDQ4kCTwJ4DgAPQF8A+ARV6PZoz2AjgDOBjAIwD8jVa8urofL1TfgwklMEXlORHqJyPkANqKld6qLDUqpIwAg8q/jH9u8RinVH8DlAG4QPScVvAzgapdjOA7AMQAWKaVWoaXVNF8pdbirUQEQkQ0i0iwiYQDPoOX8lA7WAXhbWswFEEbLhaRcp5Rqj5bW7+tux+LGKJTDIv92Q8tOeMXpGFKYAKBf5Ot+AN51MRbtKaX6oqWnWyIiW92OJ0op1SPm2ysA1LkVCwCIyGcicpiIFIpIIVqS05ki8q2bcQG7C5WoX6GlzamD8QAuAACl1PEA9oE+VyfsA6BORNa5HQhExNEHgA8B/AfAIgAXOb39mDheRctHxl1o+YO6BUAntIw+WQFgGoCOmsT1q8jXOwBsADBFk7hWAliLllbYQgBPaRLXW2hJRIsB/AvAUW7HFPfzVQDyNdlXLwL4LLKvJgA4QpO49gHwUuT/cT6AC3WIK/L88wBudzqeRA9OpSci8ijOxCQi8igmcCIij2ICJyLyKCZwIiKPYgInIvIoJnAiIo9iAici8qj/B384+c3Q98+LAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38c1ade5",
   "metadata": {},
   "source": [
    "Процедура расстановления координат по наблюдениям - Расширенный Фильтр Калмана"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d9ee774e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "m = []\n",
    "P = []\n",
    "m.append([r_t(-5),x_t(-5),y_t(-5)])\n",
    "P.append([[3,0,0],[0,3,0],[0,0,3]])\n",
    "for i in range(1,len(theta)):\n",
    "    k = theta[i-1]/T\n",
    "    m_ = np.array([m[i-1][0] + T * s_r(k), m[i-1][1] +  T * s_t(k)*cos(m[i-1][0])- T*T/2 *s_t(k) *s_r(k)*sin(m[i-1][0]),\n",
    "                   m[i-1][2] + T* s_t(k)*sin(m[i-1][0])+T*T /2 *s_t(k)*s_r(k)*cos(m[i-1][0])])\n",
    "    F_x = np.array([[1,0,-T*s_t(k) * sin (m[i-1][0]) -T*T /2 * s_t(k) *s_r(k) * cos(m[i-1][0])],\n",
    "          [0,1, T * s_t(k) *  cos(m[i-1][0]) - T*T/2 * s_t(k)*s_r(k) * sin(m[i-1][0])],\n",
    "          [0,0,1]])\n",
    "    H_x = np.array([[1,0,0],[0,1,0],[0,0,1]]);\n",
    "    p_ = np.dot(np.dot(F_x ,P[i-1]),  F_x.transpose())\n",
    "    S = np.dot(np.dot(H_x,p_),H_x.transpose())+ np.array([[sigma_r,0,0],[0,sigma_x,0],[0,0,sigma_y]])\n",
    "    K = np.dot(np.dot(p_ ,H_x.transpose()), np.linalg.inv(S))\n",
    "    v = Y[i]  - m_;\n",
    "    v = np.dot(K , v.transpose())\n",
    "    v = v.tolist();\n",
    "    m.append( m_  + v)\n",
    "    P.append(p_  - np.dot(np.dot(K ,S),K.transpose()))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "efb0323e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb710d86220>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X =[]\n",
    "for i in range(0,len(m)):\n",
    "    X.append(np.random.multivariate_normal(m[i],P[i]))\n",
    "x = [a[1] for a in X]\n",
    "y = [a[2] for a in X]\n",
    "axes1.plot(x, y,label='parametric curve')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd7cbaf8",
   "metadata": {},
   "source": [
    "Визуализация полученной траектории "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "11d4be95",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAy6UlEQVR4nO3dd3hUZfbA8e+bhCKhhiCCSAtFisIK1p8gIgoixrbqgqsYFNxV1oILWNAEG5KoiGJjlYAKrqKiiCwiKiAKCipFuqB0gYHQBkh9f3/clEmm93tnzud55kkyc+feMzfJmXfOfYvSWiOEEMJ6EqIdgBBCiMBIAhdCCIuSBC6EEBYlCVwIISxKErgQQlhUUiQPlpqaqlu2bBnJQwohhOX99NNPNq11o6r3RzSBt2zZkhUrVkTykEIIYXlKqW2u7pcSihBCWJQkcCGEsChJ4EIIYVGSwIUQwqIkgQshhEVJAhdCCIuSBC6EEBYVMwncZrORk5ODzWaLdihCCBERMZPAc3NzGTVqFLm5udEORQghIiKiIzHDKSMjo9JXIYSIdTGTwFNTUxk5cmS0wxBCiIiJmRKKMB+5LiFEeEkCF2Ej1yWECK+YKaEI85HrEkKEl7TARdiUXZdITU0Nyf6kJCNEZZLAhVdmSZxSkhGiMimhCCc2m43c3FwyMjJITU0tT5xAVHv6SElGiMqkBR4mZmm1BqJqSzcjI4Ps7OyoJ85Ql2Qixcp/C8LcpAUeJmZptQaiakvXnz72VVvvwtp/C8LcJIGHiRk+7geaTIMZFCXJypkZ/hZEbPKawJVSU4ABwD6tdefS+3KAq4ECYAuQobU+FMY4LccMI0OjkUwlWTkzw9+CiE1Ka+15A6V6AseAtx0S+BXA11rrIqXUeACt9WhvB+vevbuWVekjp6wFnp6ezuzZs6WsIYRFKaV+0lp3r3q/14uYWuvFwMEq983XWheV/rgMaBaSKEVIlbX8Zs+eLd3vhIhBoaiBDwHed/egUmoYMAygefPmITic8JcvZQ25+CiE9QTVjVAp9ShQBEx3t43WerLWurvWunujRo2COZxlRbsbmS/d72SQjBDWE3ALXCl1O8bFzcu0t0J6nLNCzwy5+CiE9QSUwJVS/YBRwCVa6+OhDSn2WCE5hqqnhJRihIgcryUUpdR7wFKgvVJqp1LqDmASUAf4Uim1Uin1epjjtDSrjiAMhJRihIgcry1wrfVAF3e/FYZYRAhEuwVshU8bQsQKmQslxkS7BRxPnzaEiDYZSh9jpAUsRPyQFniYRboLobSAhYgfksDDLNolDSFE7JISSphZoaQR7QufQojASAs8zMJZ0ghVeUY+JQhhTdICt7BQjfC0wqcEIYQzaYFbWKiWOgvFp4Roz/fiihljEiKUJIGXsuI/u5l6nJixDGPGmIQIJSmhlLLChFNmZrYyjM1mw263k5mZaZqYhAg1SeClzJaAQu0Pm53fbXYuPfPUsOzfbMuG5ebmMnbsWLKzs03xCUWIcJAEXspsCSjU3lyylXeXbeeSdo0Yc1UH2jauE+2QvAqme2OsvyELAVIDjxuPD+jEmKs68PP2PPpN/JbMT38lz14Q0RjKrjNs3LjRp+sNwdSwzXR9QIhwkRZ4rCk4DmtmQnIqnJICNepAzbpUTzqFO88s4vquPZjw1VbeWbaNT1bu5v4+bfn7BS2olhj+9/KyhLxw4ULmzp0LeL7eIK1oITzzuip9KMmq9BFwYAu8fI77x1Naw72/sPHPozw5Zx2LVm+hxtbFPPfofVx7YQe3TwvFaM2yfaSnpzN79uy4HPkpo15FIAJelV5YTP0W8MBaGPhf148f3ApfPEr7OgW8c8d5XFl9A1s+f4OMR3K4bcqPbNl/zOXTQtElr6ys0b59+4iXNyLZTdTTsaRrowglKaHEmsQkqNfMuGUdNu5blAPfPFWxzdJJsHQS6sFNPPPQvbRqVJuaHXszZcUB+r24mDsubs2/erchuUbFn4ev5YxotjCrHtvx50h2E/V0LCkLiZDSWkfs1q1bNy2iYPlbWmfWdX3b8o3Whflaa633HTmpH/xgpW4xeo4+/+kFevbKXbqkpMSvQ2VnZ2tAZ2dnBxzu/v37dXZ2tt6/f39Qx3b8OdB9BhJbqI4lRBlghXaRUyWBm0hE/vF//I/7ZJ5ZV+upV+vVq3/W/Scu1i1Gz9F/e2Op3vjnkbC8BnfbunoTcLVt1fu8/RwKoXiDEsJfksAtIFTJwWPimnq15wReeivKP67fXvqHPjvrC5328Of6yc/W6iMnCoKKqyp3r9dV/K62jUYylda1iAZJ4CESzn/gUO3ba2IrzNd6x3KtF2Z7TuTH8/SBY/n6oY9W6ZYPzdHdn/pSf/zzDr/LKu54er2+tKYlmYp4IQk8RMLR6gt1IvJrfwe2av3DZK0/Gqb1xL84J/FFOVrn2/Uv2/N0+svf6haj5+ibXv9eb957NKwxSqlCiAqSwEPE7HXVoOM7ZtN68fMuauMDdHFRsZ6+bJs+K3OebvPI5/r5+Rv1iYKigA7j7TVL6zry5Jybl7sELv3A/eTrEG1/+h07zusdbH/lQPsZlx/3hObE/t+dN/h9MQlPNmDQol4sPX8pw9NsvPHVWq6c+C3f/WbzO25vc5lHcyi8FacWDgXpo2490g88TPzpd+w4kVZOTk5Q/ZUD7WfsGO/FBVu40N2Gxw+Q/MME7gPuqwlLTpzHLW/eT/OdC/h2+osBx20m8Tq1sPRRtx5J4GES6D9DsP9Egc6q6HjcAwfSuWrECF544QXat0mDJxu6fd7FxT/yxelTmFjcmr39bqN+18spKdEkJCiPxwtnknQcwFN2LH8GFsVrIov1GTljkcyF4qd4mMui7FNAdna28Q+9dhbMvN3n5z/Z4BkGDryNNqe6nrLWZrMxadIkAIYPHx7y8+gYP1D5tcSgePibjHfu5kKRFrif4uHjtVMLtNN10LoXfDwMNs/3+vzH8h7hi0mf8E2PbG7v3cVppkN/FlsIJDm5akFHszUd7gQbD3+TwjVJ4H6KlY/XnpKKq4/SNnsxubvPI+OeV0k9vhn2rQP7Adi7Bjb+D0qKKm3fN+FHnl34Otes+zs5N55Np6b1yh/z5xwGm5zMUBYId4KNlb9JEQBXXVPCdYuFboSxwt+uix63P7TT44CgQY/k6Dc+/UafzM/36ViO3dkC6doWSLdMfwYV+Uu654lgIf3ARSCJsWy7DRs2uN++pETrLQu1/jLL6xD9n7YdDPsgnlAnfRlUJKLN0gk8Xlow4X6dgSSiQJ6Tt2qu2wS+6/GW+om7r4/oIJ5gZxC0wvQJIra5S+Bea+BKqSnAAGCf1rpz6X0pwPtAS+AP4CatdV5oizsV4uUijRlrpQHVV2fd5fahpuogjzVaQG7/YXTpfa3v+3TBsTfLwIED3a7yU3Ze7XZ7+X1Ve794qpWHs44eL3/bIkxcZXXHG9ATOAf41eG+bOCh0u8fAsZ724+WFrhXoZiKNZLcxWDbtVX/lHmRx1LKoicu1y0f+kw/NWet03B8X1v9ZdsBuk+fPhrQmZmZbuPMzMws397bVLW+vtZgmeH3KMyPYEooGC1txwS+EWhS+n0TYKMv+5EaeGB8nV410rzGUFys9W9faf3B7S6T+FMfLNItRs/RfZ5fqFfvOFT+NH/q85mZmTozM1OPHDnSbQJ3tb2/57Jsm8zMTFMkXEn88SXUCfyQw/fK8WcXzx0GrABWNG/ePHKvOIb4usBBpPkVQ942t63xTVmd9QUPv6tf/HKTLigqDn8sATy3ais+2hc0zfAGLiLHXQL3aSSmUqolMEdX1MAPaa3rOzyep7Vu4G0/sTASMxpiYaTdwc3LSZnex+3jm2t15fKDozi7WT1euKmLy1GcZjgPZojBTHGIyAj1qvR7lVJNSnfcBNgXTHDCs2jOzBcqb32ykNOeO8qapLNdPt72+Ere66vZcfA4V720hLeW/E5JSeXGRThny7PZbGRlZZGVleVxFkKz/C7CMSumsJ5AE/hsYHDp94OBT0MTTuyQf5zKMjIyeDBzPE3+MYvD177jcpvvl//KjNs60qNtKk/OWcegN5ex+9CJSvvwNAWtP+e86rZlw/vHjh0bU9OpyhSxMc5VXUVXrmG/B+wBCoGdwB1AQ+ArYDOwAEjxth8dZxcxpUbpXtm5eXvc/ZVq4RP61tDZ48frkpIS/f7y7brDY//TZ2d9oeeu3u3XfgNZY9PdBU6rM8O1EhE8gqmBh0o81cClRule2bkZ3nobp6yp3Bovrn0aic26Q53T2PaXUdz78WZW7TzM3849g8ev7kit6u6HLrg7506zK3rYVggzclcDlwQuoufkYdi2FHQxnHYWbPoC5v670iZFA17iBdv5vLZoC60aJvPSwL/Q+fR6bnZoqJqcI5Gs5Q1BhFOoL2IKEbya9aB9PzjzKqjfHNJ6O22SNOdeRv1+JzPuOJ/jBcVc9+p3TF68xekCp6Oqdd9IXHiUWrOIBkngwjwapkHWYbhvVeX7967hwoOzmHd/Dy47szHPzN3AbVN+ZO+Rky534+1ip698vShqs9mw2+1kZmbKlK4ioiSBxxHL9Ixp0NJI5I7m/pv62Y14zX4/89vPpvP2d5gwYRzf/bLG6emhanH72qou68GSnJws5RMRUbKgQxyJ1sRJ7urDXuvG962G5W/Cpnlg2wSA+nM17VjNQwmAhl2fvEP23rmMuLwdSYme2yP+1ql9nchLFlQQUeOqa0q4bvHUjTBagpkWNVxdztx17/O7q+U3zzoNw3/l7Rm6xeg5+q+vfaf3HDoRUBxCmB1uuhFKCSXGePrY76204M+FOH/KMe5q0n7VqgvssPAZp7vvzstmTrdfqLt7CYMmzmHhRveDgn05nrfX5erxaJWmLFMSE+HjKquH6yYt8MBEappZf54b8dbsrLu9rvbjeCv8Y1lAh/H2ulw9HolzYdYZKUVkEOiCDsK9SPX9nTRpEmPHjsVut5OVleVx22AWH/Dnud7qviE/N5c9Drt/NhZT9sHEWQu5JaMLp9Wr6Vds3l6XpxXvw1kDd3X9QmrvQkooQYjnvr+hLMf4pE5jGPIFNDsPqteBtMug18Nw41S44mmjT3mpbWdcS8HBndw7cTrfb3EuL/hbZnIsVQTbwyXQsoer8o9ZJtYSUeSqWR6uW6yVUCI1z4QV57OIeMwlJVqv/9yplPLumOv0/Jmv6xL7gYBjKytV9O/f3+Vz/CllhLrsYcW/DeE/rLyosYgtYU06O1a4rYsX57TV+otHtf59idaF+X7F279/f7eJN5pL4UkdPD64S+BSQhERF0x5xWsJol4zt89NOLYXvn8ZpvaHpxrBstd8Ok5qairTpk1z24PFn1JGamoqGRkZ5Obmeiyj+NobJj09vTwu6ZUSf+Qipoi4YC6+eR2MVKcxjN4GO5fD9L963lmdJoEfJwi+7NvbNq4eL5t1MRwxC3OSBB6Hoj1zni+9XdzNKJieno7dbsdut5e3jp2cUh/aXm4Mxz95hEObv6f+RzeXP/xqUTrHe4zhgQ7tSHRz/KpvMqFM6L68gZm1N4wwGVd1lXDdpAZuDlaom1aN0fFnx+99qikveMJlTXz6i6P1kUMH3D/PgVwsFNGE9AMXZazQUqsao7sWp08t495jeH/5Hqr9MpXrO1Qrv3tQ3msw4TWOtb2W2hffBc0vBKVc7iKQ/vXR/qRj9nhECLjK6uG6SQvcusLVAg12v74+32m7CZ1d91Y5ui+gOFwJ5JNOqM+z4/6s8MlLuIZ0I4xvwSaGcPVfzszM9LrfsEzC9c71rhO4/UDIpiMIZD+hPs9+l5tiRKy9VkngcS7Yf+Rw9V/OzMz0ut9A5ifxau0nWr/YpVLy/uDRAfqZ97/R454dH3ASDTYBuzrPZfdt2LAh6r83q4i1TxuSwONcOD5KR2riLG8rxgeVpDYv0HrGQF384Z2VkvlTT44N+HWFenV7x5GgsZSUwinW3rgkgYtyofrjjmQrJ9zHKnjlYtcllW3+z2oYrnJTIC1wERvcJXBZlV4ELJK9GsJ9rI/HDeX6/A/cbzByCyT7dlybzcakSZMAGD58uH+rEAnhgqxKL0IukrPheTtWsMPIew4dR07yGE62vdr1Bt++AHtWuX6sSgwAycnJjB071mm6gHiewVKEnvQDFzEh2JGS5f28Z/3D9QbLXjFupzSAdlfCFU86tcgdY8jIyHA5YtQKffCFdUgLPEZEaiKjSE6YFIpl2zztw+Vj170OD++Crn93faATebBqBsy532MMqampLlvhMoe3CCVJ4DEiUh/NI1kC8HSsqsnXXWL0tA93j9mOnuTVjQ0prtfSfXCdrnO6q2oMvq75KbMIioC5urIZrpv0QgmfWFxcwtOxfO3p4akHh7v9l+97/HitJ/f2uPZmwauXaJ1/LKjz4um1mKk7nJliiTdIN0IRSwJdVceXrn1O+y4q1Hr9HPeJ/OVz9SuZ92hAjxw50m1c7mL29Y0q2gk01gbHWIkkcBHXQpL8igrcJvHFmZdpQPfp08dtkiubNqBPnz4+x2GmuUyi/QYSz9wlcKmBi7jgWJ8OuOacWA0yD8H5zj1VerCcfVlp5D58s9e694IFC8pX5PEWh2PcvtbURfyQboQi7pRdvLTb7SQnJ/s3qEYpuHK8kcRf6lrpoUZ6Pyx+kJFdBkFKA6enDh8+vPx7n6fCdRDIlLahFM5VikRggkrgSqkHgDsBDawBMrTWJ0MRmBDhUtaCtdvtgSeklFaQdZjDJwqpN75K8l81w7j1yYKLHyi/OzU1laysLKc4rNKi9hSvjDCNjoCH0iulTgeWAB211ieUUh8Ac7XWU909R4bSCzMJVdIpXPIy1RaMcf3gP7+Hxp2iGl8klK3HmZ2dLa3zMAjXUPok4BSlVBJQC9gd5P5EDDNbf2fHkoS3uKrG7vhztYv/RdFjeZW21w3bGt+8dhG83gOWvgInD/sVXzB97iN9rqU+Hx0Bl1C01ruUUs8B24ETwHyt9fyq2ymlhgHDAJo3bx7o4USUhLIVaNYaaiCrxFf6+eqOJP13UKXt1YHNFT/8udq4Hd4F/Z7xOa5gSiyRPtfRrs/Hq4ATuFKqAXAN0Ao4BMxUSv1da/2u43Za68nAZDBKKIGHKqIh0quxR0Mgq8RX+pp4zPtBGraBXqO9blb1DTPQc27Wcx2vwlYOc9W30JcbcCPwlsPPtwGvenqO9AO3HjP1/fV3EYiIxl2YX6lf+L7s7s79xYsKvMYW7b7ekWSmv61wC/b3SqgH8gDnA2sxat8KmAb8y9NzJIGLYPjzTxC1RHhop9ZH92r9x3flifvtey4wvn9vkNaFJy0zdD7c5M3Kd+4SeFALOiilxgI3A0XAL8CdWut8d9tLL5TYY9ZFHaLeg6O4CJ5s6PKhYxePYfaK7fS99QEant46bCFE/Rx4Yfb4zMRdLxRZkUcEJV67j7lKPk73/bkGvn4KNs1zv6PHDkBieMbTxevvJhbJijwiLMzQfSwa3RNddfGret/BYydZt2UnO3pN5J+Hh7K4sEPlnfQbH/Lk7XguQvG78efcmq2baDyQofQiKGboPhaN7omuenlUve/IjKF0LPkdFt7Ha/Vc7KRVj5DHVfVcBHs+/Dm3Zu0mGsukhCIsz6y11APb1tMw9wLvGzZoBeffBZ2uhzqNgzqmv+fC2/aWuu4Qw6QGLkS05B/j2+cGcXTdAvq3rebbc1LS4JzbjFutlLCFlpWVxdixY8nMzKw0T4swF6mBi0rioV5pmte4fjY9Chf5nrwBDm6BBZmQ3Sp8cfnINOdROJEEHqciubZlpFRNNKZ5jV0GQvrL0PYK/597ifvRm6FIrMOHDyc7O7vSVLdVmeY8CidyETNOxeJQ66oX0UzzGpWqKIcAlBTD0knw5ePen1u/hduHQnHR0JeL0KY5j8KJ1MBFzLDiRbSSZ1uRcPKgy8dOdLmdU9KfM1YCcsGKr1cERmrgIuY5Lj9mFc8n/JOzX3M9Gda8DXa3yRt8e702m40p40dx+KeP4ODvxghRETOkhCJEFGUMGQJKcTxtB7VWT+OoPoU66gQA1+XPhE9rQP/noVrNgPafm5tLyYKXqHeihnFHQjVo0BIaphk9XRq2Lv2aBnWbQUJ423TyqSG0JIEL4YNwJZ7yGvSRPbB6GnXUCZ4ovJX7635D3RM74Zd3Yf0cGDwbmnTxe/8ZGRm8p49zqN951C+yGb1bDvwGB7bC1kVQdKJi48QaxlJxDdtASmuHJJ8GdZoYtfwgRXqwT6y/YUgCF8IHQS2E7Iu6TeCsG2HNTB6v9g5Lj3XkwsTSx04egjd6wqWPQvv+cOxPaNXLp2H4qamp/GtUpusHS0rg6J7SpL6l9OtWI8Fvng/FBRXbVqtlJPWyxN6wTUVyT27kc3LPyMjAbrdjt9ux2WxhT6qxPjpUErhwKdZbLv6quhByWBJ533FQlA/rZ3Nh4jrnx7952rg5anExXPca1A9gtauEBKh3unFr1bPyYyXFcHinQ3Lfanzdtw42zoUSh1p69TqVSzGOX2ulVEruqampJCcnM2rUKJKTk8OeVGO9B40kcOFSrLdc/FVW6rDZbCQnJ/u8or1fQ9FPKnK3n83dF7QjedlzbE5Mo23xFs+BbVsCL54FvR+Dnv/2+3W5lZAIDVoYt7TelR8rLoLD242EXt5y3wK7f4Z1n4Auqdi2Zr3KrfWUNIb270YN/QSDIpBUzTBXTzhJAhcuxXrLJVBVE7m38+NP6aVs24Tx43iw8VmkHdltrDbrScdrYd0n5P++lJd+yInMJ6bEpIpyStvLKz9WVMCwm/uye80Srrm4E0Ovv9RI7tuXwpqZgKY+cC/AlFznVnvDNGO/NeqE9zXECEngwiWztVzMVtLx9fxULb2A+xZ72baDMzKg/ggSlKJ4fGsSC44AUEwCiTi0bq96Hs69E4CXSuf+9rT/iEiqzoPPvM6IESPoOfoFaN++4rHCk5D3e+VW+8GtsHUhrHqv8n5qN3buJZNSmtyr14roSzIzGcgjTM9mszF48GDmzp3rtDiB2RK7O37HuW89vFoxk+GKIymMafQSn9d7lsR9v8LlT8C5Q416tH0/+ZsXMmtXKn3uHGvq8+BWgb2izl52MbUsydv3Vd627unOvWRS0oweNEk1ohN/mLkbyCMtcGF6ubm5zJ07l/79+zuVLKxSq/f7E82OHyq+f2Q3h347xqZ3VjC/Zieu5FdjGL7DUPwaQKeks0MXcKRVT4bTzjJuVZ08Uprcf6uc5NfNhhOOo1gV1DvDaLVXqbvToIXHQVFWJQlcmJ5jPb5q6zJma/VFpUvL1qgLRfn06diYsemduGzercYS4lUsqN6XfmNmMq5GbtjeyMLxacenfdasC027GrcqDuzcwtx3X+baHmdRp2BfRT/31TMh/3DFhirR6KlTtdXeMM1I+mFa1i7crBm1iCueWq++tmytUmopV1h69TL/CCx7FXqP4dame0AVO2/b/iq6Xv4S46qdG9Y3snB82gl2n1Pe+5hRD08sLa09XPGA1nD8QJV6e+nX7cugwGH6goRq0Kg9DPwv1D8j2JcUUZLARVywSqmlnGM/68U5xsXKhmmUXPwgCUueB2DDBTmcuWwkFB4ntU4NS/apDnafbp+vFCSnGrfm51d+TGs4ttdoqX/xCOxZBUUnjcFKFiMXMUVcsFwL/Oif8O3z8OPkivu6/h3OuoFjTS/ir2/8yM15k8lQnxmPqQTIzHO5K8u99kg4kQez/gGb5kGn64z52k3cdVFmIxRxLZozFQa08EKd06B/Djx2gJLqpYll5bvwznXUnnQWc0ruqUjeYAyemfcw5B8zRlBu/B8UFwKhWZAhplbl2fWzMTXBb1/BlTnw11xTJ29PpIQiRBiUtXrT09MZMWIEc+fOBQIo3yQmMS1xIBlMZnnS+Zx77rmwdJLTP25J03NIWPaqUS8vc+fX0KxbSEoflitBuaI1rJgC8x4y+pkPmQfNnBq1liItcBGTot1iLEt4ZcnbVRdIX119x2g2JHage8lPxko+tRtDn7GVtknY/bPzE798DNZ/Rmq9ZI+fPnw5VxkZGWRnZzu9hmifZ5/lH4OPh8HnI6DVJXDXYssnb5AWuIhR0W4xliW69PR0evXqFVT9OTU1ldR/vA0fZkCna+GCuyGxurHosSfbvoNt36Gr1WKdTuO0wbk0PKOt02a+nCt3vX3CPktjKOzfCO/fCgc2Q+8xcPGDYZ/3PFLkIqaISTF34e73xTDt6sr31ahXua+zB3knNJ81uIPbHppQ6X6bzcakSZMAY4Fjf89V2Xm22+2MHTvWaaRspDn93lfPhM/uM4bf3/AmtO4VtdiCISMxRVwx21wuQat7uvN9VZL3pwdacEaD6pyTsBkGTIDNX8Lxg+TXSGHK7i4MHnKn0y5yc3PLE29Z8vbnzc/fyb3CrewTQaIuYkQnGyx/E5pfaFyorNskqrGFgyRwIULM39a/T9s3TIOs0oR9eBfs+gl2rYCdP8HuX6DQzjUNt5Vvruc9jCo6CRcOp0bfp3nQzbFdLbAQSPnJLG+YGRkZ1C05xJDan8PyNXDRvXDZ4zE5jB4kgQsRcv4mQL8TZr3TsRXW4L9zfmZo45rUaNbdGIxy8lD5JqropPHN0knGvCoqwRissmeVcX+rntAni9TTuzktsGDl6QlSD6zgLv0uHNZw83ToMCDaIYWV1MBFyMVc/dlPfi3i4GcN2rHm/ErOE2we1Yz62rc6uJOUNLj356j/vkJy/OIiY7WiJS+wN+E0qt0yg5S0bqENNIpkII+ImFAMHLGisi51gM+Dhspq0MnJyT5vX9ZaH5U1HpX+ssft884cZNTDzx1afp9OrA63fgJ3LQKiO8gJQvD3cnQvvH0NLHmBVUl/ocUTm3jr469DG6RJBVVCUUrVB94EOgMaGKK1XhqCuISFWfkjeDACqR37e67Kthsy6AZqLX2Omp8+RWGtU1lZ0JruxT+iHJczAxpsmAEbZlS6b7tuSou0Sz0fSGtjdGdJMegSbLZ9TH/nbW4ZPITUxk19itVXQf29/LEEPhxiTDl77euc3uxynqyRGzd/e0GVUJRS04BvtdZvKqWqA7W01ofcbS8lFBEJ0SoJROS4WsPaWTDvIUqO7WPispN8V+0SPpoznxfGP80DrTbD2o/Dc+wyj9mif1GwpAS+nwhfPWEs7nDTO9C4Y3RjCqOQl1CUUvWAnsBbAFrrAk/JW4hIiVYJJ+yliLw/YPqN8GEGhac04j8M5PAFD/H0cy+RnZ3NrUOGwY25MGhm8Me6ZDRLqvXk0a9PsjuhogujVgl+J29/Rmv6tO2JPPjvQFiQBR2vgWELYzp5exJMCaUVsB/IVUp1AX4C7tNa2x03UkoNA4YBNG/ePIjDCeGbmCzhHNsHr14EhXZAsfPgcar/MpUrz+tF+4ONGNm3JYd3fs+Ut5ZwzaA7afjon7DwWfjuxfJd5Nc5gxpHd7je/6mdjEEuG+ZA9dpw6SN0aLOZlJLBNC1eW76ZevRPv0P3p7TkddtdP8PMwXBkjzER1XlDjalj41TAJRSlVHdgGfB/WusflFITgSNa68fcPUdKKCIcot2LIiKKCuD7l+Dg72DfR+Gh3ZywbaOOOokqKXTePiEJaqWCfR8aymvjOqGasX2PB2HvWmM6VVdung5zHnBej7J0gixfOE7oNXv2bJ975bj8XWoNK94yZlys3RhunBoTc5n4KhwjMXcCO7XWZYv3fQg8FMT+hAhItOc98SYkbzBJ1aHnv8t/rFZ6Q2ujpGDfz+Fdm1n0vw/pfV5nams72PdRcHAXf6xbTn2OcGpyQkWy//Z5z8d7/xbX908bALfMhJYXu3/u8YNw4DemvzefUaON9lxQA4LyjxlvJms+gDaXw/WToVaKT/uLdQEncK31n0qpHUqp9lrrjcBlwLrQhSaEb8xeMgnFG4zbNwGljGRWK4V6jdqT3rXywJWJOTmMeu5T+vfvz13XnEv67tK5UE7tBC0uNIaa+6PwOLx7A9z0DgeS2zLvnQmkX9iOOsd3wr51sH+DsdoNMPScYRS4mMHQ62tyFMMTUYVCsCMx/wVML+2BshUw53+QiGlmGcbtTqjn487IyPC5Re947Pr16rLz2Q9pVryDkoNbSdi31uNzXdEoY5TnjBtpCNwC8A3GcmSN2kObPlCzHix7lVqntWNk+lC3+/L6xuY4EdWtsyw7EVU4BZXAtdYrgfgpRAkRgFC8wTgmYn9a9FWPXbPv4zB3KAlFJ9DNL0Rt92/YhkKz5WAJaSlGK3hT4pmk3vI6KS27VLSM13xofG12rs+vqZKifKPWveKtmJ6IKhTks4gQFuDYRdHd4gpeFdhJbVBRO/YpebfsYZQuOl4DQN4tX3Lvlp4kP3OE3xNb0654Ayn7fqhc1ti5ApJOgcadfH5N5fK2wZS+RvK+6F4Y/BnUbWKdhSMiTBK4ECbiS6JyTHyutne7j6lXwfQbPAdww1sV3/d4EG6fAz1Hwtk3A9CgejHTxw3n86zraHzOVXDmAJg3GhbnGBdUAXatoPDUzuS88KJ/CXfjPGOtygNbjV4wVzxZ3ufcW9/+eE3wMhuhiDlW7lYYipkM3e7jjPONqWcdXK1e4cIauxhV878kNWoD1U6pePDSMXBgC+xcDt++UHrAftQHegGs/xUeWAufDoevn4L8o9DrEdizil8SujFqjI+vw2EiKk47G26aZoyudODtOoLZeyKFiyRwEXOs/M8c6Nwojtu73UfzC+GH1yvd9Zm+h92JHUk6tBWO74fN840Hap8GOa2NLoqOGnc2kuv62XDFU0YL+drXoHoyfDfRmK62uIAzr7iV7OqXeX8dR/cac5lsWwLdbod+46FaTafNvF1HMHtPpHCR6WRFzLFSCzxi86cU2OHwDnhvIOT97v05p3aE07sZFyKbnQsJifDKedD3GVg4Hpp2gdtmV4yC1NoY2l428nPEeqjrZdIrx4moBkyArgODeZUxTZZUE3HD7N0KHYXs08KJQ7D8P0aL9viB0tvBiu+L890+tZAkjqpkUsrmFR+5FZIbVt7oqDGEvuSbcegCO0cuGkMDxyHsSsHlY6H2qcYITw/J27Z/H2sn30XPooWolNZGF0EvFzyFa5LAhYgibx/9fW6h799o1KHLNDrTKHU07QK1Gla+qQSYcVP5ptUoYklRR062u4ab+l/hnLzBmB8FSCg4ylOL86mRvJiRI8933u7Cezy/4BN5HHq9P5cUb2ZDYkfOHDYfatTx/BzhliRwIaLI26cFnwfwND8f/vEdLHoW1n9mTPbU6Tq44J/GwJqqbpwKM28v/9GedhWPrGtO6x4N6O4if1M9GVQCxXVOJ7nfIG4NoNZ8aO3X8OHtpCk7C6r3pevQVyR5B0m6EQphYo59vr1Ok3taZ7j5XbjrW2jVAxaOgxfPgkXZRp3ZUcdrjS6Apa7vUJum9U7h4Y/XUFBUeVEIwCiRXPEUiTe/zQOjHvGvXq81LH+T2jP/yuG8PKZXv5U+j3xAaqNGvu9DuCQXMYWwCL8veO5ZZUwpu3Eu1KwPFw2H8+6CmnWNx4/sMS5M5h+BmvVY3P8rbnlzOd3yV/Hm0yGa1zz/GMy5H9bMpKDFJbyxvxsDh9xj+ovLZuPuIqYkcCFMKmQ9VHb/YiTyTfPglAZw0b/gvGFG+eKnafDZvcZ2Fz/ARVMOsfS9iYx+/EmeHTsmuBfgOBHVpY/IRFRBkEWNhbCYqiWTgEcbNv0LDHofhn5tdAn86gl48WxYMgE6X28MlwdYMoFpY26jcZ872FSvO+4adz7FsXomTL4UThw0epn0HCnJOwzkjAphUlXnPAl0qbjyhFujhTGX951fwennGP22J3aB1Lbl27bdMo0Xnnqcn/cV8+FPO13uz2McRfkwZwR8fCc0Oduox8ssgmEjvVCEMKmqPVQCHW3o1Ne8WXf4+0ew40fjQueKKRUb/zSVv919D7NaNuDpuevpfeapNKxdo9L+3MaRt81Y7mz3L8ZEVJc9Hv3Fj2Oc1MCFsJBA6uJen7P9B/j6Sfjj2/K7tgzbzD0fbGDCzV3p0KSu94NsnAez7jJ6nFz7KnQY4P05wmdSAxciBgRSRnFsybusXTc/35h18KqKZdbSJrflf+f/SodUzy1o274/Wfb0lfDezVC/Ody1UJJ3BEkCF8JCAp4LHPfJv7xG3upauODu8vvVFw/DxK7wwxtQeNJ5h0f3cvz1y7mg8HtWJf0F7vjSaRZBEV5SAxfCQoKZ58Vd7bpSjfyBLFj2qvHAecPgz1/hf6NgyYvQYwSccxsk1SifiOqMhMPMTUrnvGETXM4iKMJLauBCRIlZZk10iuOPJcbiDwCP58Efi+GbcbBjGdRpasxMeHgHNGwDN70tE1FFgNTAhQiRUK3+Emi3wFApex1A5aXNWl5sDPgBY7Wd1r1gyDxjbcqju43kDXDOYGjY1nnHImIkgQvhp1Al3mDq2Z74+gbj8XXct9r4+uNk4+vuX2BBJiRUg3b94PTu8OVj8HI3+GkqFBeG9DUI30gNXAg/hWr1l3DNW+7rHOMeX0fNujDoA2M+leVvGqvE125stMSbdTe6C/62AL55Bj67D7593hht2WWg9P2OIKmBC2FiYen37auiAvj0HljzAbS5HK6fDLVSKm+jtbEM28JxRiu9QcvSRZD/BonSPgwVqYELYUHB9PsO+sLo/g1G8gbQJcZkWCcOVd5GKWjXF4Z+AwP/a8w9/uk9MKk7rJxhLFgswkZa4EKYWNR7quxeCWs/hrWz4NB2owbe5jLodD20v7JiatoyWhvT1y4cB3+ugZQ0uGQUnHWj0XtFBESmkxVCBE5r2PVzaTL/BI7shMQa0KaPMaNhu76VV9fRGjbMMaax3fur0eXwktHQ+QZJ5AGQBC5EDItoS72kBHatMFrlaz8xuhYm1YS2lxst83Z9jSXYyrbd8JmRyPetg9R2RiLvdJ0kcj9IAhcihuXk5DBq1Ciys7PD0rPFrZIS2PGDkczXfQLH9kLSKUYS73y9cfGzei1ju/WfwsLxsH+9sejyJaONpd1knnCvJIELEWGRbBVHvVYOUFIM25eWJvNPwb4fqiVD+35Gy7xNH0isDutmGYncthFO7Wgk8g7pksg9kAQuRIRlZWUxduxYMjMzycrKinY4kVVcBNu+M5L5+tlw/ABUr2Nc+Ox8vTG6c8PnRmnlwGZo3NlI5GcOkETugnQjFEJETmIStL4Ern4RHtxkLKvW+Tr47Ut472/wXDtjINAVT8I1r0LRSfjgVnijJ6yfY1wEdSNUUxnEAulpL0SYDB8+nOTk5JAPlbecxCRI623crnoBti4yerOsnwOr3jP6jrfvD03PMSbMev8WOO1s6PWw0WJXqtLufB1pGg+CLqEopRKBFcAurbXHmdylhCJEZaaoXUdLUT5s+cYos2z4HAqOQo16UJxvtMgBmnQ1Enm7vuWJPB7PWThLKPcB60OwHyHiTrRnJIyqpBrYGnYnZ0sbbEN+gL/NgHZXQIJDYWDPSmO1n//0hk3zQevQjTSNAUGVUJRSzYCrgKeBESGJSIg4EqqJsazKqRxy5lVQeAI2f2m0zDfNg8LjsPtnmHGjMQtir4eN0aBVSivxKKgSilLqQ2AcUAf4t5RQhBD+8FoOKTgOm78oTebzoeiEcX+zc6H3GKM3SxwIeTdCpdQAoL/W+m6lVC/cJHCl1DBgGEDz5s27bdu2LaDjCSHiXP4xo0W+dpbRQi8phJFbnGdIjEHhSODjgFuBIqAmUBf4WGv9d3fPkRa4ECIkTh6BvN+hSZdoRxIRIb+IqbV+WGvdTGvdEvgb8LWn5C2EECFTs27cJG9PZCCPEEJYVEgG8mitFwILQ7EvIYQQvpEWuBBCWJQkcCGEsChJ4EIIYVGSwIUQwqIkgQshhEVJAhdCCIuSBC6EEBYlCVwIISxKErgQQliUJHAhhLAoSeBCCGFRksCFEMKiJIELIZzYbDZycnKw2WzRDkV4IAlcCOEkrhdbtpCQTCcrhIgt8b7YslVIAhdCOElNTTVWiRemJiUUIYSwKEngQghhUZLAhRDCoiSBCyGERUkCF0IIi5IELoQQFiUJXAghLEpprSN3MKX2A9uAVMCMY3QlLv9IXP4xY1xmjAkkrqpaaK0bVb0zogm8/KBKrdBad4/4gb2QuPwjcfnHjHGZMSaQuHwlJRQhhLAoSeBCCGFR0Urgk6N0XG8kLv9IXP4xY1xmjAkkLp9EpQYuhBAieFJCEUIIi5IELoQQFhXxBK6Uuk8p9atSaq1S6v5IH98hjilKqX1KqV8d7ktRSn2plNpc+rWBSeK6sfR8lSilotKFyU1cOUqpDUqp1UqpWUqp+iaJ68nSmFYqpeYrpZpGOyaHxx5USmmlVGokY3IXl1IqSym1q/RcrVRK9TdDXKX3/6v072utUirbDHEppd53OFd/KKVWRjouRxFN4EqpzsBQ4DygCzBAKdUmkjE4mAr0q3LfQ8BXWuu2wFelP0faVJzj+hW4Hlgc8WgqTMU5ri+Bzlrrs4FNwMORDgrXceVorc/WWncF5gCPmyAmlFJnAFcA2yMcT5mpuIgLmKC17lp6mxvhmMBFXEqpS4FrgC5a607Ac2aIS2t9c9m5Aj4CPo5CXOUi3QLvAPygtT6utS4CFmEkpojTWi8GDla5+xpgWun304BrIxkTuI5La71ea70x0rFUicFVXPNLf48Ay4BmJonriMOPyUBEr9S7+dsCmACMinQ8ZTzEFVVu4von8KzWOr90m30miQsApZQCbgLei2hQVUQ6gf8K9FBKNVRK1QL6A2dEOAZPGmut95R+/yfQOJrBWMwQ4H/RDqKMUupppdQO4BYi3wJ3Fc81wC6t9apox+LC8NKS05RolA3daIeRK35QSi1SSp0b7YCq6AHs1VpvjmYQEU3gWuv1wHhgPjAPWAkURzIGX2mjf6X0sfSBUupRoAiYHu1YymitH9Van4ER0/BoxlLaWHkEE7yRuPAakAZ0BfYAz0c1mgpJQApwATAS+KC01WsWA4ly6xuicBFTa/2W1rqb1ronkIdROzWLvUqpJgClXyP+sc1qlFK3AwOAW7Q5BxVMB26IcgxpQCtglVLqD4xS089KqdOiGhWgtd6rtS7WWpcA/8G4PmUGO4GPteFHoARjIqmoU0olYZR+3492LNHohXJq6dfmGCdhRqRj8GA2MLj0+8HAp1GMxfSUUv0warrpWuvj0Y6njFKqrcOP1wAbohULgNZ6jdb6VK11S611S4zkdI7W+s9oxgXlDZUy12GUOc3gE+BSAKVUO6A65pmdsA+wQWu9M9qBoLWO6A34FlgHrAIui/TxHeJ4D+MjYyHGP9QdQEOM3iebgQVAikniuq70+3xgL/CFSeL6DdiBUQpbCbxukrg+wkhEq4HPgNOjHVOVx/8AUk1yrt4B1pSeq9lAE5PEVR14t/T3+DPQ2wxxld4/FfhHpONxdZOh9EIIYVEyElMIISxKErgQQliUJHAhhLAoSeBCCGFRksCFEMKiJIELIYRFSQIXQgiL+n86w2SFyzxpgQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88af2c8e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "001cd767",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
